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ABSTRACT 

Classical Be stars are an enigmatic subclass of rapidly rotating hot stars characterized by dense equatorial 
disks of gas that have been inferred to orbit with Keplerian velocities. Although these disks seem to be ejected 
from the star and not accreted, there is substantial observational evidence to show that the stars rotate more 
slowly than required for centrifugally driven mass loss. This paper develops an idea (proposed originally by 
Hiroyasu Ando and colleagues) that nonradial stellar pulsations inject enough angular momentum into the up- 
per atmosphere to spin up a Keplerian disk. The pulsations themselves are evanescent in the stellar photosphere, 
but they may be unstable to the generation of resonant oscillations at the acoustic cutoff frequency. A detailed 
theory of the conversion from pulsations to resonant waves does not yet exist for realistic hot-star atmospheres, 
so the current models depend on a parameterized approximation for the efficiency of wave excitation. Once res- 
onant waves have been formed, however, they grow in amplitude with increasing height, steepen into shocks, 
and exert radial and azimuthal Reynolds stresses on the mean fluid. Using reasonable assumptions for the 
stellar parameters, these processes were found to naturally create the inner boundary conditions required for 
dense Keplerian disks, even when the underlying B-star photosphere is rotating as slowly as 60% of its critical 
rotation speed. Because there is evidence for long-term changes in Be-star pulsational properties, this model 
may also account for the long-term variability of Be stars, including transitions between normal. Be, and shell 
phases. 

Subject headings: circumstellar matter — stars: early-type — stars: oscillations (including pulsations) — stars: 
emission-hne, Be — stars: rotation — waves 



1. INTRODUCTION 

Be stars are non-supergiant B-type stars that exhibit emis- 
sion in their hydrogen B aimer lines. It has been known for 
many decades (e.g., Struve 1931) that the "classical" Be stars 
tend to rotate more rapidly than normal non-emission B stars. 
A wide range of observations of Be stars is consistent with 
the coexistence of a dense circumstellar disk in the equatorial 
plane and a variable stellar wind at higher latitudes (Doazan 
1982; Slettebak 1988; Prinja 1989; Porter & Rivinius 2003). 
There is a great deal of evidence that the so-called decretion 
disk is ejected from the star and is not accreted from an exter- 
nal source of gas. The physical mechanisms that are respon- 
sible for producing the disk, though, are still not known. This 
paper investigates the idea that nonradial pulsations (NRP) 
may deposit sufficient angular momentum above the photo- 
spheres of rapidly rotating Be stars to explain the origins of 
dense equatorial disks. 

Theoretical explanations for the Be phenomenon depend 
crucially on how rapidly the stars and disks are rotating. Al- 
though there is increasing evidence that the extended disks 
are in Keplerian orbit (e.g., Hanuschik 1996; Hummel & 
Vrancken 2000; Okazaki 2007), there is also evidence to show 
that the underlying photospheres rotate too slowly to propel 
any atmospheric gas into orbit. Typical observational values 
of the ratio of equatorial rotation speed to the critical rota- 
tion speed Vciit (at which gravity is balanced by the outward 
centrifugal force) range between 0.5 and 0.8 (Slettebak 1982; 
Porter 1996; Yudin 2001). There has been some recent dis- 
agreement between this traditional picture and new ideas — 
inspired by the tendency for gravity darkening^ to mask rapid 
rotation at the equator — that most or all Be stars may be ro- 

' An oblate, rigidly rotating star is expected to undergo an internal re- 
distribution of its net radiative flux in proportion to the magnitude of the 
centrifugally-modified gravity, which produces hotter and brighter poles and 



tating very nearly at V^ni (Townsend et al. 2004). Cranmer 
(2005) found that this may be true for the late-type Be stars 
(i.e., spectral types B3 and later), but early-type Be stars do 
seem to be consistent with a range of intrinsic rotation speeds 
between ^40% and 100% critical. For stars rotating suffi- 
ciently below Vciit, any physical model for the origin of Be- 
star disks requires a significant increase in angular momentum 
above the photosphere. 

This paper develops a set of ideas concerning how stellar 
NRP may give rise to outwardly propagating circumstellar 
waves, which in turn can deposit angular momentum in Be- 
star disks. The history of this theoretical perspective is dis- 
cussed in more detail in § 2. However, a number of other disk 
formation mechanisms have been proposed, and it is impor- 
tant to keep in mind that in some cases these could be acting in 
cooperation or competition with the pulsational processes ad- 
vocated here. Some of the other suggested mechanisms have 
been: (1) stellar "wind compression" that may channel super- 
sonically outflowing gas from mid-latitudes to the equatorial 
plane by conservation of angular momentum (Bjorkman & 
Cassinelli 1993; Bjorkman 2000), (2) episodic ejections from 
some point on the star, with some material being propelled 
forward into orbit and some propelled backward to fall onto 
the star (Kroll & Hanuschik 1997; Owocki & Cranmer 2002), 
and (3) magnetic forces that can channel gas toward the equa- 
torial plane and provide sufficient torque to spin up a disk 
(e.g., Poe & Friend 1986; CassinelH et al. 2002; Brown et al. 
2004, 2008; Ud-Doula et al. 2008). 

The three mechanisms listed above have one general fea- 
ture in common: they depend on the existence of large-scale 
supersonic flows in the circumstellar regions. Thus, these 
processes appear likely to give rise to strong variability on 

a cooler, dimmer equator (see, e.g., von Zeipel 1924; Slettebak 1949; Collins 
1963; Maeder 1999). 
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short (dynamical) time scales. Although many Be stars do ex- 
hibit variability on time scales of hours to days, many other 
Be stars seem to maintain their disks in a relatively quies- 
cent manner over several years (see, e.g., Dachs 1987; Telting 
2000; Okazaki 2007). It is worthwhile, then, to first consider a 
"baseline" disk formation mechanism that feeds the disk with 
a more or less steady-state supply of angular momentum from 
below. Other proposed mechanisms then may act as sources 
for the complex observed patterns of Be-star variability (see 
also § 6.2). 

The outline of this paper is as follows. § 2 presents an 
overview of how stellar pulsations may be coupled to the out- 
ward transport of mass and momentum in the circumstellar 
regions of Be stars. In § 3 the detailed properties of waves 
and shocks are presented, including how evanescent pulsation 
modes couple to propagating waves, how these waves steepen 
into dissipative shocks, and how the waves can exert a net 
pressure to affect the time-steady properties of the stellar at- 
mosphere. § 4 describes the time-independent conservation 
equations for mass and momentum that are solved with contri- 
butions from the waves, and § 5 gives the results. A discussion 
of the implications of these models on our overall understand- 
ing of Be-star disk formation and time variability is given in 
§ 6. Finally, § 7 contains a brief summary of the major results 
and discusses how the proposed physical processes should be 
further tested and refined. 

2. THE DISK-PULSATION CONNECTION 

Do Stellar pulsations represent a causal factor in producing 
the wide variety of observed manifestations of the Be phe- 
nomenon? This possible connection has been discussed for 
several decades both from the standpoints of Be-star observa- 
tions (e.g., Baade 1982, 1985; Henrichs 1984; Willson 1986; 
Smith 1988) and pulsation theory (Ando 1986; Osaki 1986; 
Castor 1986; Saio 1994; Lee 2006, 2007). For a while, there 
was significant debate over whether much of the observed 
Be-star line profile variability was even the result of NRP, or 
if these features were simply rotationally modulated "spots" 
(e.g., Baade & Balona 1994). However, recent improvements 
in the detection and modeling of NRP in many Be stars has 
put the pulsational interpretation on firmer footing (Rivinius 
et al. 2003; Rivinius 2007; Townsend 2007). 

The well-studied Be star /i Cen has been a key target in 
the search for disk-pulsation correlations (e.g., Rivinius et 
al. 1998). The multiple NRP periods of /x Cen seem to un- 
dergo beating in phase with outbursts of circumstellar mate- 
rial into the disk, suggesting a direct input of mass and mo- 
mentum when the pulsation displacements are largest. In- 
creases in disk emission at NRP amplitude maxima have also 
been observed for other multiperiodic Be stars such as ( Oph 
(Kambe et al. 1993a) and 28 Cyg (Tubbesing et al. 2000). 
However, similar kinds of outbursts are also seen for Be stars 
Uke oj CMa that seem to have only a single pulsation period 
(Stefl et al. 2003). For other Be stars, there has been line 
profile variability observed at Doppler velocities that exceed 
the projected photospheric rotation velocity (Chen et al. 1989; 
Kambe et al. 1993b). This may be evidence for a radial in- 
crease in the rotation rate between the photosphere and the 
inner edge of the disk. 

In addition to observational connections between pulsations 
and disks, there are many stars without disks (O, B, and Wolf- 
Rayet) for which there exists coupled variability between the 
photosphere and the stellar wind. The radial pulsator BW 
Vul clearly shows nonlinear shock-like features in its pho- 



tosphere (e.g.. Smith & Jeffery 2003) that are in phase with 
radially accelerating features in its supersonic wind (Massa 
1994; Owocki & Cranmer 2002). Several stars show different 
periods for photospheric and wind variabiUty, with the former 
(due to NRP) often being shorter than the latter (e.g., Howarth 
et al. 1993, 1998; Kaufer et al. 2007). Although much of the 
observed stellar wind variability may be attributed to struc- 
tures corotating with the star (Kaper et al. 1996), the dis- 
tribution of these stars in luminosity and effective temper- 
ature seems to agree roughly with the domain of so-called 
"strange-mode" oscillations (FuUerton et al. 1996), indicating 
that there may be a pulsational cause for some types of wind 
variability. 

From a theoretical standpoint, it has been known for 
some time that small-scale hydrodynamic fluctuations (waves, 
shocks, turbulence, pulsations) can affect the large-scale mean 
properties of a fluid. Linear waves propagating in an inho- 
mogeneous medium exert a time-steady wave pressure on the 
fluid even without dissipating (Bretherton & Garrett 1968; 
Dewar 1970; Jacques 1977). Alfven wave pressure is be- 
lieved to be a key contributor to the acceleration of the high- 
speed solar wind (e.g., Cranmer 2004). A nonlinear exten- 
sion of this effect is likely to be acting in the outer atmo- 
spheres of cool pulsating giants and supergiants — e.g., Mira 
variables and asymptotic giant branch stars — with pulsation 
driven mass loss suspected to occur in many cases (Willson 
2000; Woitke 2007; Neilson & Lester 2008). 

In addition to supplying linear momentum, it has been 
shown that waves can transport and deposit angular momen- 
tum in rotating systems. Interior gravity waves in rotating 
planetary atmospheres have been shown to be responsible 
for driving various kinds of mean shear flows (Eliassen & 
Palm 1960; Bretherton 1969; Plumb 1977; Schatzman 1993; 
Rogers et al. 2008). In stellar interiors, the study of how diffu- 
sive transport processes interact with the mean rotation profile 
has proceeded with the impUcit assumption that waves or tur- 
bulent eddies may be the cause of the "anomalous" transport 
(e.g., Rudiger 1977; Zahn et al. 1997; Talon 2008). Ando 
(1982) considered the coupling between waves and rotation 
as a possible means of generating differential rotation in stel- 
lar interiors. Recently, Townsend & MacDonald (2008) found 
that the interior transport of angular momentum by pulsations 
may be an important factor in the evolution of rotating mas- 
sive stars. 

The apphcation of the above ideas to the episodic spinup 
and mass loss of Be stars was developed by Ando (1986, 
1988), Osaki (1986, 1999), Lee & Saio (1993), Saio (1994), 
and Lee (2006, 2007). This mechanism was originally be- 
lieved to transport angular momentum up to the stellar sur- 
face only for prograde NRP modes (i.e., modes for which 
the azimuthal group velocity is in the same direction as the 
mean rotation); see, e.g., Ando (1983). A retrograde oscilla- 
tion mode was thought to transport angular momentum down- 
ward and thus potentially spin down the outer layers of the 
star. However, more recent work — which takes a more com- 
plete inventory of the various second order wave-rotation cou- 
pling terms — finds that some types of coupling do not depend 
at all on the sign of the azimuthal group velocity (e.g., Lee 
2007). Furthermore, Townsend (2005) and Pantillon et al. 
(2007) showed that some "mixed modes" may exist in rapidly 
rotating stars that appear to have retrograde phase velocities 
but prograde group velocities. There are thus several potential 
ways for retrograde NRP modes (which appear to be preva- 
lent in Be stars) to be associated with the transport of angular 
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momentum up to the surface and beyond. 

Aside from the work of Ando and colleagues, though, 
surprisingly little work has been done to study how pulsa- 
tional energy (as well as the effects of pulsations on the 
mean fluid) may "leak" out of a hot star into its circum- 
stellar medium. Cranmer (1996) studied some aspects of 
how low-frequency evanescent waves in B-star photospheres 
might evolve into propagating waves at larger heights in a stel- 
lar wind. Townsend (2000a,b,c) modeled NRP leakage from 
the perspective of whether the p and g modes are propagating 
or evanescent in the layers directly beneath the photosphere 
(see also Gautschy 1992). The goal of this paper is to follow 
up on these ideas and to investigate how an extended hot-star 
atmosphere responds to the presence of pulsational oscilla- 
tions at its photospheric base. 

3. WAVES AND SHOCKS 

This section contains a derivation of some of the physical 
processes that determine how pulsationally driven waves can 
influence the mean properties of a stellar atmosphere. The 
subsections below isolate four of the main effects that need 
to be considered: evanescence (§ 3.1), resonant excitation at 
the acoustic cutoff (§ 3.2), shock steepening and dissipation 
(§ 3.3), and wave pressure (§ 3.4). Some of the proposed con- 
nections between the above mechanisms are still somewhat 
speculative, so it is hoped that this combination of ideas will 
eventually be tested with numerical simulations. 

This paper makes the implicit assumption that there is a rea- 
sonably time-steady input stream of pulsational oscillations 
from the stellar interior. Observations suggest, however, that 
many NRP modes have finite lifetimes, and that they undergo 
secular amplitude variations over a wide range of time scales 
(see § 6.2 for discussion of the implications of NRP transience 
on Be-star variability). The models presented below also ig- 
nore several effects that may be important in some situations, 
such as heat conduction, radiative damping, magnetic fields, 
and ballistic freefall behind strong shocks (the latter being im- 
portant in, e.g., Mira variables). 

3.1. Evanescent Pulsation Modes 

For the purposes of modeling the effects of B-star NRP os- 
cillations on their outer atmospheres, it is important to first 
summarize what is known about them observationally. The 
f3 Cep variables have typical pulsation periods of order 2-10 
hrs, spectral types between BO and B2, and are often inter- 
preted as low-order p and g modes (Sterken & Jerzykiewicz 
1992; Aerts & de Cat 2003). Slowly pulsating B (SPB) stars 
have longer periods of 10-50 hrs, later spectral types of B3 
to B9, and they appear to be pulsating in high-order g modes 
(e.g., de Cat 2007). Pulsating Be stars often have periods that 
appear similar to the SPB class, but with some modes (includ- 
ing beat periods) that extend up into the 100-200 hr range 
(e.g., Kaufer et al. 2006). Some B stars are "hybrid" pulsators 
that exhibit both (3 Cep and SPB type modes (Dziembowski 
& Pamyatnykh 2008). 

Figure 1 shows a collection of B-star pulsation frequencies 
oj and horizontal wavenumbers that have been normalized, 
respectively, by dividing by an estimated acoustic cutoff fre- 
quency uja in the photosphere and multiplying by an estimated 
density scale height H. Spectral types, periods, and angular 
mode identifications were obtained for the /3 Cep and SPB 
stars from de Cat (2002, 2003)^ and for Be stars from Riv- 
inius et al. (2003). Additional low-frequency modes for two 



1.0 r 

I A 




0.002 0.004 0.006 0.008 0.010 



Fig. 1 . — Collected B-star pulsation frequencies and azimuthal wavenum- 
bers, normalized to dimensionless quantities via the photospheric acoustic 
cutoff frequencies Wa and scale heights H. Data for /3 Cep variables (red 
triangles), SPB stars (green squares), and Be stars (blue circles) are shown. 
Filled symbols assume nonrotating stellar properties; open symbols attempt 
to take account of rapid rotation. Also shown are acoustic-gravity wave prop- 
agation boundaries (solid lines) and /-mode curve (dotted line). 

hybrid B-type pulsators (12 Lac and ly Eri) were taken from 
Dziembowski & Pamyatnykh (2008), and the assignment of 
£=\m\=2 to these modes is relatively uncertain. 

The two-dimensional iuj,kx) plot illustrated by Figure 1 
is often called a diagnostic diagram for waves in stellar at- 
mospheres. The solid curves are the "propagation boundary 
curves" between which acoustic-gravity waves are evanes- 
cent. Above the upper curve, waves propagate vertically as 
acoustic or p-mode waves; below the lower curve, they prop- 
agate as internal gravity or g-mode waves (e.g.. Lamb 1908, 
1932; Mihalas & Mihalas 1984). To determine the normaliz- 
ing values of the acoustic cutoff frequency coa and scale height 
H for each star, the tabulated spectral types and luminosity 
classes were converted into stellar masses M», radii R„, and 
effective temperatures Tgff in a single uniform way; i.e., using 
the luminosity and Teff relations of de Jager & Nieuwenhui- 
jzen (1987), and interpolating onto the evolutionary tracks of 
Claret (2004) to obtain masses (see also § 2 of Cranmer 2005). 
The magnitude of the horizontal wavenumber was computed 
from the tabulated meridional degree number £ and the as- 
sumption that ( = \m\, the latter being the azimuthal mode 
number. Thus, \k^\ = \m\/R^,. For the Be stars, only the modes 
identified as having m = 2 by Rivinius et al. (2003) are shown. 

The solid symbols in Figure 1 show values computed using 
the standard expressions for uja and H for nonrotating stars 
(see below). The open symbols show modified values that 
take account of two attempts at improving the accuracy of 
these normaUzations: 

1 . The surface gravity was reduced by taking account of 
both centrifugal and radiative acceleration. A lower 
limit to the stellar rotation was appUed by assuming that 



The 2005 version of this database was obtained from: 



|http://www.ster.kiileuven.ac.be/^peter/Bstars/| 
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the observed projected rotation rate vsin/ is equal to the 
actual equatorial rotation speed Vrot- For the Be stars, it 
was assumed that the rotation speed is never smaller 
than Vcrit/2 (see Cranmer 2005). The radiative acceler- 
ation in the photosphere was estimated using a fit to the 
numerical results of Lanz & Hubeny (2007) given by 
equations (l3TI)-(l32b below. 

2. For the Be stars only, the observed pulsation frequency 
was replaced by an estimate of the intrinsic pulsation 
frequency in the rotating frame, assuming that the Be- 
star oscillations are retrograde (Rivinius et al. 2003). In 
this case, the rotating-frame frequency lu is higher than 
the observed inertial-frame frequency cjq, with 

oj LO{) + m^{\—C„i) (1) 

and C„i « 0.2 for low-order / and g modes (e.g., 
Ledoux 1951). We ignored second order centrifugal 
terms in the above expression (i.e., terms proportional 
to Vl^ /bj) since for low-order g modes with I = \m\ = 2 
they would give a negligibly small correction (Saio 
1981). 

For the SPB stars, these corrections were found to be unim- 
portant and are not shown. 

Taking account of the above effects results in smaller val- 
ues of LOa and in larger values of both lo and H, thus moving 
each data point upwards and to the right in Figure 1 . Even 
with these corrections, though, nearly all of the observed pul- 
sation modes are still solidly evanescent (i.e., not propagating 
vertically) in the photosphere. Also, the Be-star data points 
seem to have moved from the upper edge of the SPB region 
into the (3 Cep region. (Many of them seem to cluster around 
the fundamental-mode curve, given by uP' = gk^, which is a 
rough boundary between the /9-mode and ^-mode regions.) 
This may suggest that classical Be-star pulsations have a sim- 
ilar intrinsic character as those of the jS Cep variables, despite 
the apparent similarity of their periods to those of the SPB 
stars (see also Baade 1982). 

The observed ranges of oscillation parameters can be put 
into context by using a theoretical model of acoustic-gravity 
waves. The model presented here makes the simplifying as- 
sumptions of an isothermal (constant temperature T) medium, 
plane-parallel geometry, and a constant gravitational acceler- 
ation g. It is also straightforward to first consider the system 
in the rotating frame, in which we define the frequency ui. In 
this atmosphere, the sound speed is constant and is given by 



(2) 



where kg is Boltzmann's constant, nin is the mass of a hy- 
drogen atom, and /i is the mean molecular weight of the gas. 
The latter quantity is assumed to be /i = 0.604, which corre- 
sponds to a standard (10% helium) elemental abundance mix- 
ture. The adiabatic exponent 7 is the ratio of specific heats 
Cp/cy, and usually can range between 1 (isothermal) and 5/3 
(adiabatic). For an isothermal atmosphere, the pressure and 
density scale heights are equal to one another and given by 
H = c^/'yg (which does not depend on the value of 7 itself). 

The standard assumption in linear wave theory is to model 
the physical variables (i.e., density, velocity, pressure) as the 
sum of a time-steady (zeroth order) component and an oscil- 
lating (first order) component. The linear amplitudes vary as 
^iut-tk,z-ihx ^jjjj ^ frequency lo. We model the conditions 




in the equatorial plane, with gravity acting in the -z direction 
and the stellar rotation proceeding in the +x direction (see Ap- 
pendix). The azimuthal wavenumber can thus be expressed in 
conventional NRP terminology as k^ = -m/Rf,, with negative 
[positive] values of m corresponding to prograde [retrograde] 
modes. The time-steady vertical velocity wq is assumed to be 
zero, and the azimuthal velocity mq is also set to zero (in the 
rotating frame). The zeroth-order density po is proportional to 
g-z/H ^ The corresponding linear amplitudes are denoted wi. 
Ml, and p\, and they are modeled as complex quantities in or- 
der to account for phase differences between their respective 
oscillations. 

For a given frequency u) and horizontal wavenumber k^, the 
vertical wavenumber is in general a complex quantity given 
by 



2H 



■± 



(3) 



where the acoustic cutoff frequency is Ua = ^g/2cs and the 
Brunt-Vaisala frequency is ujg = (7 - 1)'/^^/^ (see, e.g.. 
Lamb 1932; Mihalas & Mihalas 1984; Wang et al. 1995). The 
regions of {ijj,k^) space where the vertical wavenumber has a 
real part correspond to propagating waves, and the regions 
where fe, is completely imaginary correspond to evanescence. 
If the amplitude and phase of the vertical velocity amplitude 
wi are specified, the horizontal velocity and density fluctua- 
tions are given in dimensionless form by 



Ml 
Wi 



ctk. 



k-- 



Pi/ A) 

Wl/Cs 



-cW. 



UJC, I I 

k, + — 
H 



{l-l)c]kl 



■1 



(4) 



(5) 



Note that when k^ is completely imaginary, the above expres- 
sions are also imaginary. This implies that in the case of 
evanescent waves, mi and w\ are 90° out of phase, and p\ 
and w\ are also 90° out of phase. 

It is useful to illustrate how these waves behave in the limit 
of small horizontal wavenumbers (i.e., \k^H\ -C 1) which ap- 
plies to B-star NRP observed at the stellar surface (see Figure 
1). It is also useful to use the simplifying limit of 7 = 1, which 
corresponds to isothermal fluctuations. This approximation is 
consistent with the existence of rapid radiative heating and 
cooling in the relatively dense outer atmospheres of hot O 
and B stars (e.g.. Drew 1989; Millar & Marlborough 1999; 
Carciofi & Bjorkman 2006). With the above assumptions, the 
lower propagation boundary disappears and the upper propa- 
gation boundary is greatly simplified: i.e., waves with ijj> Ug 
propagate; waves with oj <LOa are evanescent. 

Propagating waves have exponentially growing amplitudes 
with height, with the non-oscillating part of the height depen- 
dence being given by w\ oc e'l^'^ . These waves have a vertical 
phase speed 



ph 



Re A:, 



1-4 



-1/2 



(6) 



and a group velocity Vgr = c^/Vph. For upwardly propagat- 
ing waves, the above amplitude relations can be simplified to 
show that 



Re(i^ 



2k,H^^ «1 
UJ Vph 



(7) 
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(8) 



Note that the sign of k^- determines the sign of the magnitude 
ratio Mi/wi, but the ratio pi/wi is positive definite. 

Evanescent waves have two possible solutions for 
their (completely non-oscillatory) height dependence; one 
"steeper" than the propagating solution, and one "shallower" 
The shallow [steep] solution has a total wave energy density 
that decays [grows] with increasing height. Thus, the shal- 
low solution is the physically realistic choice when consid- 
ering waves that originate at low heights and have effects on 
the medium at larger heights (see, e.g., Wang et al. 1995). The 
vertical amplitude dependence of the shallow solution is given 
by 



wi oc exp 




(9) 



For evanescent waves, the phase speed is formally infinite 
and the group velocity is zero. However, Cranmer (1996) 
showed that if a time-steady subsonic wind (wo ^ c?) is in- 
cluded in the linear equations for acoustic-gravity waves, the 
shallow evanescent solution can be shown to have a large — 
but finite — upward phase speed of order c^/wq. The steep 
evanescent solution corresponds to a downward phase speed 
of -c^/3w(). These both correspond to negligibly small (sub- 
sonic) group velocities and do not need to be considered fur- 
ther. 

3.2. Resonant Wave Excitation 

It has been known for some time that the acoustic cutoff uja 
is a preferred resonance frequency for an atmosphere that has 
been disturbed by a pulse or piston-like initial condition (see, 
e.g.. Lamb 1908; Schmidt & Zirker 1963). The acoustic cut- 
off frequency also acts as a fundamental "ringing" mode for 
disturbances of a sinusoidal nature (Fleck & Schmitz 1991; 
Kalkofen et al. 1994; Sutmann et al. 1998). Of particular 
interest here is the case of an evanescent oscillation driven 
by NRP (with uj < cj,,) that can excite a resonant oscillation 

(W « UJa).^ 

The one-dimensional problem of resonant excitation in 
an infinite isothermal atmosphere can be solved analytically 
(e.g.. Fleck & Schmitz 1991). In this section we make use of 
this simple model to gain insight about the physics of resonant 
waves. Although the model contains features that are formally 
inapplicable to actual pulsating stars (i.e., the assumption of 
an instantaneous "start" at f = 0), it remains valuable as a way 
to obtain numerical estimates of the resonant wave properties 
in the absence of a more accurate theory. Let us then consider 
a piston with an average height slightly below the photosphere 
(Zpiston < 0) with a vertical pulsation amplitude Vpuis. For sim- 
plicity, the long-wavelength horizontal velocity fluctuations 
can be ignored; i.e., ks = 0. Thus, for a piston oscillating with 
evanescent frequency uj, the height dependence of this driven 
mode is given by equation (|9]l, and its time dependence at 
larger heights remains strictly sinusoidal. In addition, when 
the piston motion begins (at time f = 0) there arises a resonant 

^ In addition to long-period evanescent waves giving rise to resonant exci- 
tation, it is important to note also that shorter period acoustic waves may also 
lead to resonant excitation at the cutoff, possibly via nonlinear effects such 
as "shock cannibalization" (see Rammacher & Ulmschneider 1992). Thus, 
there appear to be several independent ways to create waves at the cutoff 
frequency. 



response at the acoustic cutoff frequency. These two compo- 
nents of the vertical velocity amplitude can be expressed as 



Vi-es(z,f) 



sm Loj 



wi(z,f) = 

T 



T/ -iujt+nZ , 



H COS 



3tt 

iOat- — 



(10) 



where Z = {z-Zp\ston)/'^H and k = l-(l-cj^/[x)^)'''^ (Kalkofen 
et al. 1994; Sutmann et al. 1998). In the late-time limit of 
t ^z/cs, after which the initial transient disturbance has prop- 
agated away, the resonant component's velocity amplitude is 
given by 



Vres(z,t) = Vpuls 



Ze^ 



TT (a;2-w2)f3/2 



(11) 



This expression gives only the lowest order term in an infinite 
expansion; the higher order terms have more rapid time decay 
(e.g., f"^/^, f"''/^, etc.) and thus the above term dominates the 
late-time (f z/c,) behavior of the resonant mode. 

For a given time f after the piston begins to oscillate, both 
the driven (evanescent) mode and the resonant (cutoff) mode 
coexist in the atmosphere. Immediately above the piston, 
Vres Vpuls, but as one proceeds higher in the atmosphere, 
the radial growth of the resonant mode overtakes that of the 
shallow evanescent mode (eq. ||9l) and the resonant mode be- 
comes stronger. It is suspected that the stellar photosphere sits 
below the height where this transition occurs, because the ob- 
servable signatures of NRP show only the evanescent driver 
A major conjecture of this paper is that the upper layers of 
the atmosphere — including the "inner edge" of the Keplerian 
disk — are those where the resonant oscillations are dominant. 

Equation (fTTT l shows that the resonant oscillation amplitude 
decays in time as f"^/^. If left undisturbed it would decay 
eventually to zero, leaving only the driven evanescent mode. 
However, this property of the solution seems to be an artifact 
of the idealized nature of an infinitely extended, isothermal, 
and plane-parallel atmosphere. For a real atmosphere, there 
are several effects that have been suggested to arrest this decay 
and produce a quasi-steady stream of resonant waves: 

1. Fleck & Schmitz (1991) found that a sharp reflecting 
boundary at some height above the piston (even a rela- 
tively permeable one) can cause enough downward re- 
flection of the resonant waves to maintain them at a fi- 
nite amplitude over long times. Reflection effectively 
"restarts the clock" on the time t in equation ( fTTT l every 
time the piston sees a reflected wave. 

2. Radial gradients in the sound speed or scale height 
give rise to gradual linear reflection that can lead to 
continual excitation of the resonant waves (e.g., Pitte- 
way & Hines 1965; Schmitz & Fleck 1992; Lou 1995; 
Musielak et al. 2005; Erdelyi et al. 2007; Taroyan & 
Erdelyi 2008). This kind of reflection is also believed 
to be acting on Alfven waves in the solar chromosphere 
and corona, and in that case it may be responsible for 
seeding an ongoing turbulent cascade (Heinemann & 
Olbert 1980; Matthaeus et al. 1999; Cranmer & van 
Baflegooijen 2005). 

3. If the driven evanescent oscillations are multiperiodic 
or stochastic in nature, the resulting departures from 
simple sinusoidal piston motion could result in con- 
tinuously excited resonant waves. Such intermittency 
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could be the result of rotationally generated subsurface 
convection (Espinosa Lara & Rieutord 2007; Maeder et 
al. 2008), NRP mode beating (Rivinius et al. 1998), or 
other nonlinear processes that can cause one NRP mode 
to grow at the expense of another mode (see, e.g.. Smith 
1986). 

The remainder of this paper assumes the existence of a con- 
tinual regeneration of resonant oscillations in B-star atmo- 
spheres. The general phenomenon of wave reflection seems 
to be the most likely route for this to occur. Fleck & Schmitz 
(1991) showed how the presence of wave reflection may cause 
the resonant waves to "forget" about an instantaneous initial 
condition and become essentially self-sustaining. Although 
the reflection in the models of Fleck & Schmitz (1991) was 
induced by a sharp upper boundary, there are several other 
ways that reflection can be produced gradually in an actual 
stellar atmosphere (see the second item in the above list). In 
the absence of a robust theory or simulation of these effects in 
hot-star atmospheres, however, we are limited to using sim- 
pler prescriptions such as equation ( fTTT l. 

The use of the simple theory derived above should not be 
interpreted as a suggestion that resonant waves are produced 
by an abrupt initial condition. This model is used only to 
estimate the mean amplitude of a steady-state ensemble of 
resonant waves, by choosing a representative time t in equa- 
tion ( fTTT l. Equivalently, this can be specified in terms of the 
height A of a fictional reflecting boundary above the photo- 
sphere, and the reflection time is thus estimated as f « A/c, 
(see, e.g.. Figure 2 of Fleck & Schmitz 1991). The param- 
eter A is used as a free parameter in the models discussed 
below. Eventually, of course, there should be time-dependent 
numerical simulations of the partial wave reflection and reso- 
nant wave excitation in B-star atmospheres that would allow 
the most appropriate values of A to be determined. 

What values of A should be considered realistic? It is use- 
ful to estimate the heights in the atmosphere at which the fluid 
parameters are thought to undergo qualitative changes. These 
are the heights at which substantial wave reflection should oc- 
cur One suggestion is to compute the height at which the 
strongest wave modes steepen into shocks (and probably sat- 
urate in amplitude). This onset of nonlinearity signals a rapid 
change in the overall properties of the atmosphere. For hot 
stars with NRP amplitudes that are akeady of the same order 
as the sound speed in the photosphere, the height at which 
they steepen into shocks may be as low as 2 to 4 photospheric 
scale heights above the photosphere itself. For a star with 
an accelerating wind, other significant heights would be the 
sonic and critical points (which are not the same for a radia- 
tively driven wind; see Castor et al. 1975). For the standard 
B2 V star model discussed in detail in § 5, the sonic point 
would be at a radius of about 10 scale heights. The critical 
point is another scale height or two above the sonic point. 

It is possible, though, that the relevant time scale f could be 
much larger than suggested above. For a B star, a value of 
A « \QH corresponds to a time that is only about 10% of a 
single evanescent NRP period. If t is determined instead by 
long-term secular changes in the NRP "piston" (i.e., beating 
or mode growth/decay), it may be as large as the mode life- 
times themselves; typically no smaller than 10 to 100 NRP 
periods. However, since there is also evidence for shorter time 
scales being important (consistent with ongoing wave reflec- 
tion from the upper atmosphere), we will make the order-of- 
magnitude assumption that A is of order 10 scale heights and 




z / H 



Fig. 2. — Comparison of height dependences for driven evanescent wave 
amphtudes Vp^is (dashed line), resonantly excited wave amphtudes Vrcs (solid 
lines), and the maximum available wave amphtudes Vmax from energy conser- 
vation (dotted lines). Thick lines denote ideal linear results; thin hnes show 
numerical results that take account of wave dissipation and shock steepening. 

vary it up and down from that baseline value. 

Figure 2 compares the height dependences of vertical veloc- 
ity amplitudes for evanescent and resonant waves. The stellar 
parameters assumed here are those of the baseline B2 V model 
discussed in detail in § 5. The evanescent NRP driver has a pe- 
riod of 10 hrs, the acoustic cutoff period is 1 .3 hrs, and the pis- 
ton is assumed to sit 0.5// below the photosphere. For an as- 
sumed value of A = 10//, the resonant mode becomes stronger 
than the evanescent mode at a height of approximately 4 scale 
heights above the photosphere (i.e., when the hydrostatic den- 
sity has dropped by a factor of ~50 from the photospheric 
value). Also shown is an upper-limit velocity amplitude V„,ax 
that illustrates the total energy density available for waves in a 
stratified atmosphere. This quantity was computed under the 
assumption that poVj^ax constant, and V^ax is normalized by 
the evanescent NRP driving amplitude at the piston. The thick 
curves show linear undamped solutions (i.e., eq. ifTTI ) and the 
thin curves show the result of applying the shock steepening, 
dissipation, and wave pressure effects discussed below. 

Note that the radial growth of the undamped resonant mode 
is more rapid than for a classical acoustic-gravity wave at 
w > (the latter illustrated by Vmax fx e^). The extra lin- 
ear factor of Z in equation (fTTl i is a unique feature of the 
resonantly excited oscillation. This factor has the strongest 
impact very near the piston height, but grows progressively 
less important once one goes a few scale heights above the 
piston. The resulting properties of the upper atmosphere are 
thus relatively insensitive to the exact value chosen for Zpiston- 
Also, when wave damping and wave pressure are taken into 
account (i.e., the thin curves in Fig. 2) the resonant mode no 
longer grows faster than V^ax at all. 

Choosing larger values of the reflection height A results in a 
lower "efficiency" in the conversion of driving wave energy to 
resonant wave energy. This efficiency is estimated as the ratio 
Vres/Vmax, and Figure 3 shows how this varies as a function of 
the A parameter itself. The models that were used to compute 
these efficiencies contained the full range of shock steepen- 
ing and wave pressure effects (see §§ 3.3-3.4), and they con- 
strained Vies to never exceed the maximum available energy 
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Fig. 3. — Efficiency of resonant excitation (Vres/Vmax) plotted as a function 
of the assumed reflection height A for a series of B2 V atmosphere models 
(see text). The photospheric {dotted line) and maximum (solid line) values 
for each model are shown. 

specified by Vmax- The efficiencies are shown both at the pho- 
tospheric base and at the height of peak efficiency. The latter 
occurs usually a few scale heights above the photosphere, be- 
fore substantial damping sets in. It is clear that values of A/// 
less than about 3 or 4 may not be physically realistic, since 
they result in the resonant wave attempting to extract more 
energy from the driving mode than is actually available. 

At the subphotospheric height of the idealized piston, the 
models below ignore the effects of radiative acceleration giad 
(which can counteract gravity at and above the photosphere). 
Lanz & Hubeny (2007) reported the existence of a local max- 
imum in grad at the r w 1 photosphere, with less of an effect 
in the deeper layers that have smaller temperature gradients. 
Thus, because the net gravitational acceleration is higher at 
the piston than in the photosphere, the acoustic cutoff fre- 
quency at the piston is assumed to be slightly higher than that 
at the photosphere. For the standard B2 V star model, this 
results in a resonantly excited wave having a frequency about 
1 .2 times the local acoustic cutoff frequency once it reaches 
the photosphere. In other words, the assumption is that the 
wave eventually "forgets" about its origin as a resonant ex- 
citation and just propagates (as any other linear wave would) 
with a constant frequency. 

3.3. Wave Action Conservation and Shock Dissipation 

The resonantly excited wave modes discussed above (uj > 
LOa) are assumed to propagate upward from the photosphere 
and grow in amplitude. At some point, these mainly longi- 
tudinal fluctuations are expected to steepen into shocks and 
begin to dissipate some of their energy (e.g.. Castor 1986). 
The nonlinear evolution of such a wave train is modeled here 
using a formalism described more completely by Cranmer et 
al. (2007). The wave energy density Us is assumed to obey an 
equation of wave action conservation. 



d_ 
di 



]_d_ 



(wo + Vgi)At/.5 



(12) 



(see also Jacques 1977; Koninx 1992), where A is the cross- 
sectional area of a flow tube and 2sh is the shock dissipation 



rate. The time derivative term in equation ( fT2] i is given only 
for completeness and is ignored in the time-stationary models 
below. In general, the Doppler shifted frequency in the mov- 
ing frame (i.e., taking account of both rotation and a stellar 
wind) is given hy lo' = uj- uok^ - wqL. As above, the mq effect 
is ignored by interpreting cu as the corotating frequency, and 
the factor dependent on wo is taken into account by using the 
proper dispersion relation. The wave energy density is defined 
as 

Us = -pw\ (13) 



where s is a dimensionless shape factor determined by the 
spatial profile of the waves. The peak-to-peak vertical veloc- 
ity amplitude is denoted w\ , and the variance w\ is equivalent 
to the magnitude of the product w\w\ (see Appendix). 

The adopted area function A(z) contains both the spherical 
expansion of the stellar atmosphere and an additional linear 
factor (depending on the distance above the piston) that de- 
scribes how the resonantly excited waves differ from standard 
acoustic-gravity waves (eq. IfTTI ). Thus, 



A (X 



(z- 



" ^piston )^ 



(14) 



where r is measured from the center of the star, and thus only 
varies by a few percent over the heights considered here. The 
models use the assumption that Zpiston is 0.5// below the pho- 
tosphere. Although A diverges unphysically at z = Zpiston, this 
location is always kept outside the modeled grid of heights. 
The precise normalization of A is unimportant, and only its 
relative variation with height enters into equations ( fT2l i and 
O- 

The Doppler shifted frequency oj' can be computed in terms 
of the comoving-frame phase speed Vph, and a more useful 
version of the wave action conservation equation can be writ- 
ten as 



d_ 
dr 



(W0 + Vgr)(W0+Vph)At/, 



(wo + Vph)A6sh 
Vph 



(15) 



At the photospheric base, the lower boundary condition on the 
wave energy density is specified as Us = with s = 2 

(i.e., a sinusoidal wave shape) also assumed at this height. 

Before proceeding to describe the shock steepening in more 
detail, it is useful to examine the wave action equation in a 
simpler limiting case. For a negligibly small vertical wind 
speed Wo, and for waves sufficiently above the acoustic cutoff 
frequency that Vph « Vgr ~ Cj, equation (fTSl l simplifies greatly 
to become 

c.?^ = -esh, (16) 

or 

where the additional assumption of A w constant was also 
applied for large heights z ^ Zpiston in a plane-parallel at- 
mosphere. This equation has the benefit of clearly illustrat- 
ing how shock dissipation affects the radial gradient of the 
wave energy density. The other terms that have been stripped 
away are non-dissipative contributors to the wave energy gra- 
dient. Although the numerical solutions presented in § 5 uti- 
lize the more complete equation of wave action conservation 
(eq. ifTSll ). the simplified form of equation ( fTSI l is later found 
to be useful to derive a physicaUy meaningful expression for 
the wave pressure. 

A standard way of expressing the time-averaged shock dis- 
sipation rate is 

_pT^S 



(17) 
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where AS is the net entropy jump across a shock, and it 
is nonzero only above the height where the wave train has 
steepened into shocks. This expression uses an approxima- 
tion from so-called "weak-shock theory" that the volumetric 
heating rate is given by the internal energy dissipated at one 
shock divided by the mean time between shock passages in a 
periodic train. This assumption breaks down for very strong 
radiative shock trains (e.g., Carlsson & Stein 1992, 1997), 
which dissipate their energy in relatively narrow zones behind 
each shock. However, the models presented below do not de- 
velop such strong shocks. The gain in internal energy across 
an ideal inviscid shock is given by 



TAS 



c,, 



72-7^1 (Pi/Pi 



(18) 



where c,. is the specific heat at constant volume for an ideal 
gas, and subscripts 1 and 2 denote quantities measured on the 
upstream (supersonic) and downstream (subsonic) sides of the 
shock (Landau & Lifshitz 1959). The above expression is 
inconvenient for the limit of isothermal shocks (7=1) since 
the term in square brackets approaches zero and the specific 
heat Cy is proportional to (7-1)"' and thus diverges to infinity. 
Their product remains finite, and can be rewritten in this case 
as 



TAS = ct 



M\-\ 
2M? 



-InM^ 



(19) 



Az 



(20) 



where Mi is the Mach number of the shock. The above ex- 
pressions are valid for shocks of arbitrary strength and are not 
limited to the traditional weak-shock approximation. 

In order to use the shock dissipation rate in the wave ac- 
tion conservation equation, the Mach number M\ needs to be 
expressed in terms of the local velocity amplitude w\. Cran- 
mer et al. (2007) described a method of following the steep- 
ening of an initially sinusoidal wave profile {s = 2) into a 
fully steepened sawtooth or N-wave (s = 3). Those steepen- 
ing equations are modified here slightly due to the fact that 
Vph y Cj for acoustic-gravity waves having vertical wavelength 
A = iTrV^h/u). The instantaneous distance between a wave 
crest and the zero-velocity node immediately ahead of it is 
given by 

Ao 7+1 /" w\dz 
T 2~y 'V^ ' 

and the integration is taken from the photospheric lower 
boundary (at which A = Ao) up to a given height. The progres- 
sive nonlinear steepening of a wave gives rise to a decrease 
in Az as the wave train propagates upwards. At every point 
above the photosphere, the local steepening factor C = Az/A 
is known and is used to compute both the shape factor s and a 
shock efficiency e = (Mi - l)cj/wi. The latter quantity is the 
ratio of the shock velocity amplitude to the full velocity am- 
plitude of the wave profile, and it can range between zero (for 
a sinusoidal wave that has not yet steepened) and one (once 
the crest overtakes the node ahead of it and C, < 0). Cranmer 
et al. (2007) gave both numerical results and analytic fitting 
formulae for siQ and e(Q. The models presented in § 5 make 
use of these analytic fits. 

3.4. Wave Pressure 

Waves that propagate up from the stellar photosphere, 
steepen into shocks, and dissipate a fraction of their energy 
also are able to transfer a mean momentum flux (i.e., exert 
a ponderomotive force) to the bulk atmosphere. Although in 



some situations it is possible for dissipationless waves to ex- 
ert pressure (or Reynolds stresses) in a fluid, the ability for 
dissipating waves to do so is much easier to understand (e.g., 
Goldreich & Nicholson 1989; Koninx 1992). The oscillatory 
energy lost via dissipation is sufficiently "randomized" to pro- 
duce a net, time-averaged source of both linear and angular 
momentum. 

In the Appendix, the second-order wave pressure contri- 
butions to the radial and azimuthal momentum conservation 
equations are computed from first principles. This derivation 
does not make use of the standard linearizing assumption that 
the first-order oscillations are smaller in magnitude than the 
corresponding time-steady (zeroth order) equilibrium proper- 
ties. It does, however, assume that the variations of the first- 
order quantities in time (?) and in the azimuthal direction (0 or 
x) are sinusoidal. T he res ulting momentum source terms are 
given in equations ( lAllI ) and dAOI l. and the wave-pressure 
acceleration is thus expressed as 



-(S)/po 



(21) 



where the vector g„p is defined as positive on the right-hand 
side of the radial momentum conservation equation (eq. 1301 ) 
when it points away from the star 

The radial component of the w ave-pressure acceleration is 
given by applying equation ( IA13b . with 



1 dUs 
^"P''' " po dr 



(22) 



Note that Us is constant in the limit of small-amplitude 
acoustic-gravity waves that propagate (i.e., with uj > LOg) in 
a hydrostatic, isothermal, and plane-parallel atmosphere. In 
that case, gwp.r = 0. Dissipation is thus required for a net up- 
ward wave-pressure acceleration, and the wave action conser- 
vation equation derived in § 3.3 gives the required value of 
dUs/ dr. Equation ( fTSI l can be used to obtain a good approxi- 
mation for the dissipation rate, and 



>TAS TAS 



Sv/p.r — 



4ttH 



(23) 



The use of equation ( fTSI l implies that only the direct shock 
dissipation is being considered as a contributor to a nonzero 
wave energy gradient. It may not be correct to apply the more 
complete equation of wave action conservation (eq. |[T5l ) to 
the expression for wave pressure, since the additional radial 
derivatives of A, Vgi, and Vph are not dissipative in origin. 
These other contributors to dUs/dr are essentially there to 
keep the wave energy flux locally conserved, and they do not 
represent irreversible energy exchange with the surrounding 
medium. 

Note that the right-most approximation in equation ( |23T l 
makes the assumption that u k, lu^, which is motivated by 
the fact that propagating waves are suggested to arise from 
resonant excitation (see § 3.2). The above expression for 
gv/p.r thus depends inversely on the density scale height. In 
a subsonic atmosphere, one of the effects of wave pressure 
is to produce a radial increase in H as the wave amplitude 
increases. As // ^ oo, the background medium would be- 
come homogeneous, and it would make sense for gwp.r ^ 
because a nonzero wave pressure gradient relies on there be- 
ing a background inhomogeneity. Thus, it is not clear whether 
the scale height to be used above should be the value at the pis- 
ton (which is consistent with the assumption of constant fre- 
quency w), or if the locally varying value H(z) would be more 
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consistent with the known behavior of wave pressure. Be- 
cause the proposed Be-star spinup mechanism depends cru- 
cially on the magnitude of gwp, we aim to be safe and thus un- 
derestimate the local value of gwp.r by using the larger value 
of the local scale height H{z) in the denominator of equation 

das. 

The azimuthal component of gwp is not often discussed in 
wave-pressure acceleration literature (e.g., as applied to stel- 
lar winds), but it is a key component in studies of the interac- 
tions between waves and mean flows in planetary atmospheres 
(see § 2). Equation (lAl lb gives 

1 



^wp,0 



' d 

— (Po(miWi)) + MOTr((PlM'l)) + 

po [or Or 



{PlWl 



duf) 
dr 



(24) 



The first term, which depends on the phase-averaged prod- 
uct {u\W\), takes the form of a Reynolds stress that resembles 
the radial wave pressure (i.e., with a factor of ui replacing 
one of the two factors of w\ implicit in eq. 1221 ). The sec- 
ond and third terms depend on the phase-averaged product 
{p\W\), and they represent a radial transport of "eddy mass 
flux" (Lee 2006). In the theoretical NRP literature, the {uywy ) 
and {p\Wi) effects are associated with angular momentum 
transport time scales tj and To respectively (e.g., Lee & Saio 
1993). 

Ando (1983, 1986) found that the ti transport term gives 
a net outward transport of angular momentum for prograde 
modes (m < 0, or kx > 0) and inward transport for retro- 
grade modes. This can be understood from the perspective 
of acoustic-gravity wave theory by examining the sign depen- 
dence of equation (|7]l- The tj transport term can only spin up 
the outer layers of a star when m < 0. The T2 transport term, 
however, does not depend on the sign of m at all, and may 
also be nonzero for purely radial pulsations (m = 0) in a rotat- 
ing star. Lee (2006, 2007) found that for rapidly rotating stars, 
the T2 term is much stronger than the ri term near the stellar 
surface. This can be verified by using the approximate scal- 
ing relations given in § 3.1 to estimate the relative strengths 
of the two effects. The ratio of the first of the two T2 terms in 
equation (|24] | to the ti term is given roughly by 



PlUQ 

PoUl 



UJUq 
C^kx 



1 



> 1 



(25) 



2k, 

The third term in equation (|24] | acts somewhat like an effec- 
tive viscosity. It was derived under the assumption of a plane- 
parallel geometry to be proportional to du^ / dr, but in a spher- 
ical geometry it should in fact be proportional to the full radial 
gradient of the angular momentum ruQ. This generalization is 
included in equation (|26] |. 

In the models presented below, the amplitude and phase re- 
lations derived in § 3.1 are used to scale the azimuthal wave 
pressure with the radial wave pressure. Thus, the expression 
used is 



,^ ^ X {piwi) f duQ Mo 

gv.p.4> = (*i + ^'2)^wp,,-^^^ -[ 1^ + — 

Po \dr r 



where 



$1 = 



(uiwi) 



Mo(piWl) 



Poiw]) 



(26) 



(27) 



Note that when acoustic-gravity waves are evanescent, both u i 
and wi are 90° out of phase with one another (making <i>i = 0) 



and p[ and W[ are 90° out of phase (making $2 = 0). When 
applying these transport terms to stellar interiors, Ando (1983, 
1986), Lee & Saio (1993), and Lee (2006, 2007) depended 
on nonadiabatic effects to produce departures from these 90° 
phase differences. In the present case, however, these subtle 
effects do not need to be invoked. Radially propagating waves 
in stellar atmospheres naturally exhibit $1 t^^O and $2 7^0. 

4. CONSERVATION EQUATIONS 

In order to evaluate how the wave interactions derived in § 3 
impact the overall atmosphere, both the mean and fluctuating 
fluid properties must be computed simultaneously. The time- 
averaged density and flow speeds obey conservation equa- 
tions for mass, radial momentum, and angular momentum 
that are given below. These equations are solved for time- 
independent fluid properties with the implicit assumption that 
the wave pressure terms — despite having their origin in tem- 
porally fluctuating motions — are essentially time-steady. The 
energy conservation equation is assumed to be satisfied by a 
known and constant temperature T. 

The mass conservation equation for a spherically symmet- 
ric atmosphere can be written in two ways, depending on 
whether one tracks the mass flux through a fixed location in 
the star's inertial frame, or whether individual fluid parcels 
are followed in time. The former (Eulerian) treatment yields 
the most commonly seen version of the equation of mass con- 
servation. 



M = Airpowor 



(28) 



where M is the mass loss rate averaged over the entire star 
and r is measured from the center of the star. Alternately, 
equation ( IA7b can be used to determine the more complete 
(Lagrangian) mass flux, with 



M = 47r(/9ovvo+(piM'i))r^ 



(29) 



The proper physical interpretation of equation ( [29] ). however, 
is not clear In the photospheres of pulsating stars (including 
the Sun), the magnitude of the so-called "Stokes drift" term 
(piWi) may exceed the time-steady mass flux M/(47rr^) by 
at least an order of magnitude. For example, in the unified 
photospheric, chromospheric, and coronal models of Cran- 
mer et al. (2007), the maximum of the Stokes drift velocity 
(pi wi ) /po for the Sun occurs in the upper chromosphere with 
a magnitude of ^7 km s"' (about half the local sound speed), 
and at this height the mean upflow speed of the solar wind 
(computed effectively by solving eq. 1281 for wo) is only 0.4 
km s"'. If equation ( |29] ) were the correct time-steady mass 
conservation equation, it would imply that the Eulerian radial 
velocity wq would have to be negative in these regions — i.e., 
approximately equal in magnitude to the Stokes drift velocity, 
but opposite in sign — to nearly cancel out the Stokes drift and 
produce the smaller known mass flux. This appears to be a 
kind of unphysical "fine tuning" that tends to appear in situ- 
ations where some key piece of physics has been neglected. 
On the other hand, equation (|29] | does seem to be a more self- 
consistent treatment of the second-order wave pressure terms 
as derived in the Appendix. 

In the models below, we take an agnostic approach to this 
dilemma and simply choose the version of the mass flux con- 
servation equation that produces a smaller angular momen- 
tum transport in the upper atmosphere of a Be star. If this 
choice was the wrong one to make on the basis of physical 
realism, then correcting it would only increase the ability of 
NRP/wave coupling to be an effective agent in spinning up 
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Be-star disks. The more "conservative" choice was found to 
be equation (|28] ) rather than equation ( l29] l: see below. The 
former equation is solved with the mass loss rate M speci- 
fied with a known constant value and the density stratification 
poif) constrained by the radial momentum equation. Equation 
(|28] ) is thus solved for the local value of the Eulerian outflow 
speed wq. 

The radial momentum equation is given in the subsonic 
Umit (i.e., wq <^ c j as 



7Po dr 



GM^ ui 
— — + - + 



(30) 



where the first two terms on the right-hand side (gravity and 
the centrifugal force) are given in their full spherical form. As 
discussed above, we limit ourselves to isothermal fluctuations 
having 7=1. The radiative acceleration gr^d is specified using 
the numerical model atmosphere results of Lanz & Hubeny 
(2007). The photospheric value of g,ad is assumed to apply 
for the entire atmospheric model, and thus the additional de- 
pendences on the density and the velocity gradient that apply 
to stellar winds (Castor et al. 1975) are not applied here. 

Figure 4 illustrates a parameterized fit to the Lanz & 
Hubeny (2007) radiative acceleration results as a function of 
stellar gravity and effective temperature (compare with their 
Figure 19). The so-called Eddington-limit acceleration ^edd is 
the surface gravity at which a star of a given T^ff has sufficient 
radiative acceleration to completely cancel out its gravity. The 
effective temperature dependence of this quantity has been fit 
with a power law, with 



^edd = 142.35 



'eff 



2 X 104 K 



(31) 



The resulting radiative acceleration, expressed as a dimen- 
sionless ratio to the surface gravity (g > 0), was then fit with 
the following polynomial relation. 



In 



^rad 



-0.25943ji:-0.037106x'' 



(32) 



where x = ln(g/§edd)- This expression is valid only for g > 
^edd (i-e., below the Eddington limit). Figure 4 also shows 
how \ogg depends on Tgff for standard definitions of stars 
having luminosity classes 111, IV, and V, using the spectral 
type calibrations of de Jager & Nieuwenhuijzen (1987) and 
Cranmer (2005). In the Be-star models below, we use the 
centrifugally-modified gravity (i.e., the first two terms on the 
right-hand side of eq. ||30l ) for g in the parameterization for 

gad- 
In practice, it is simpler to solve equation ( l30b for the local 
scale height H than to integrate it directly to obtain the density 
Pq. Using equation ( |23] l for the wave-pressure acceleration, 
the scale height is given by 



H 



c^ + (TAS/4tt) 



where the effective gravity 



^eff = 



CM, 



^rad 

r 



(33) 



(34) 



is always positive. Note that the scale height can become 
large in two interesting limiting cases: (1) when the shock 
dissipation TAS is substantial, and (2) when the atmosphere 
spins up enough so that the effective gravity geft is small. 
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Fig. 4. — Contours of constant values for the ratio of radiative to gravita- 
tional acceleration, as a function of stellar surface gravity g and effective tem- 
perature Teff, computed from fits to the numerical models of Lanz & Hubeny 
(2007). Values of the ratio giad/g are listed next to each curve, with the Ed- 
dington limit (g^-iilg = 1) shown as a thicker curve. Also shown are locations 
of luminosity classes III (dotted line), IV (dot-dashed line), and V (dashed 
line) stars in g, T^fi space. 

The azimuthal component of the momentum conservation 
equation, which determines the bulk rotation speed mq, is 



duQ M() 



(35) 



No horizontal radiative acceleration is assumed here, though 
it may be present in the supersonic regions of rapidly rotat- 
ing winds or disks (e.g., Grinin 1978; Owocki et al. 1996; 
Gayley et al. 2001; Kubat 2007). Equation ^ is used for 
the azimuthal wave-pressure acceleration gv,p.ci>, and thus the 
conservation equation can be rearranged as 



(poWo+{piw 



. / duo Mo\ _ 



($i-l-$2)pogwp.r . (36) 



It is clear that without any horizontal wave pressure, the equi- 
librium solution for the rotation speed would be mq oc r"' ; i.e., 
angular momentum conservation and no spinup to a Keple- 
rian disk. The mass flux given in the first set of parenthe- 
ses above is identical to the Lagrangian terms from equation 
(IA7I) . When deciding which form of the mass flux conserva- 
tion should be applied in these models (i.e., eq. 1281 or 1291 ) 
it was concluded that the safer choice would be to use the 
version that gives the larger mass flux. In that case, the right- 
hand side of equation ( |36] ) would be divided by a larger quan- 
tity, and the radial gradient of the angular momentum would 
be increased by a smaller (i.e., more conservatively underes- 
timated) amount. As mentioned above, this choice was found 
to be equation ( |28] |. since in that case both wq and the Stokes 
drift velocity {p\W\) / pq are always positive. 

Equations ( [15] ), ( l28T l. ( l33b . and ( |36] | were solved simulta- 
neously by a numerical code that integrates upwards from 
the photospheric lower boundary. The numerical quadrature 
was done with straightforward explicit Euler integration steps. 
The density po was advanced along the radial grid by using 
the local value of the scale height H. Most runs of the code 
utilized 4000 radial grid zones that covered 80 photospheric 
scale heights, and tests with double the number of grid zones 
verified that the resolution was adequate. 
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Because the radial momentum conservation (eq. 1301 ) does 
not account self-consistently for the transonic and supersonic 
parts of the stellar wind, additional checks are made in the 
code to ensure that the subsonic assumption is not violated. If 
the integrated value of wq exceeds a specified terminal speed 
Woo, then the wind speed is not allowed to increase any fur- 
ther and the density scale height is set to an asymptotic value 
^^oo = r/2. These solutions are judged to be "winds" and not 
"disks." Furthermore, the wave amplitude wi is monitored 
so that it does not exceed V^ax (which scales with and 
is normalized to the photospheric boundary condition for wi). 
This ensures that the modeled resonant excitation does not ex- 
tract more energy from the evanescent waves than is available. 

5. MODEL RESULTS 
5.1. Example B2 V Star at 70% Critical Rotation 

Figure 5 shows the result of integrating the conservation 
equations derived above for an example B2 V star rotating 
well below its critical rotation speed. For this star, M» = 
7.43 M0 and logL*/L0 = 3.47, and the polar radius was as- 
sumed to be Rp = A.HRq (see Cranmer 2005). The equa- 
torial rotation speed Viot was chosen to be 70% of the crit- 
ical speed, with Vait = 479 km s"' for this star. Assuming 
rigid rotation and Roche equipotentials, the equatorial radius 
is Rf = A.92Rq. Standard von Zeipel (1924) gravity darken- 
ing gives the equatorial Tgff to be 17920 K, which is lower 
than both the corresponding polar value of 22170 K and the 
nonrotating value of 20900 K. The equatorial scale height at 
the photosphere is Hq = 0.00177?^, and the maximum height 
shown in Figure 5 is 8O//0, or ~0.14 R^. 

The photospheric boundary condition on density po is com- 
puted from a tabulated grid of Rosseland mean opacities kr 
(Kurucz 1992) and the condition that r w ky^pqHq = 1. The 
value for the example star, at the equator, was computed to 
be approximately 7 x 10~'° g cm"^. This model uses a repre- 
sentative B2 V star mass-loss rate M = 10"^ M0 yr"', which 
assumes that the polar wind's mass-loss rate (estimated here 
using the fitting relations of Vink et al. 2000) remains com- 
parable to the equator's mass loss rate (see § 5.1 of Bjork- 
man 2000). Also, we assume a stellar wind terminal speed 
Woo = 1500 km s~'. Note, however, that the adopted values 
of M and Woo tend not to have any real impact the quantita- 
tive properties of "successful" disk models — i.e., models that 
do not exceed the supersonic wind condition (wq = Woo) — 
because the wq term in equation (|36] | is dominated in those 
cases by the viscosity-like (Stokes drift) term. 

The example stellar model shown in Figure 5 was assumed 
to have a driving NRP period of 10 hours, which is equiva- 
lent to an evanescent frequency ratio oj/tJo = 0.13 in the pho- 
tosphere (see Fig. 1). This frequency enters into the model 
calculations only in the factor of (uol-Lip-) in the denominator 
of equation ( fTTT i. so the model is relatively insensitive to the 
exact value of the driving frequency. Elsewhere in the model, 
the wave frequency lo is taken to be that of the resonantly ex- 
cited waves at the subphotospheric piston. The acoustic cutoff 
period at the piston is 1 .07 hours, and this is slightly shorter 
than the local acoustic cutoff period in the photosphere (the 
latter using grad to lower the effective gravity), which is 1.32 
hours. The adopted photospheric NRP amplitude is 20 km s"' 
(i.e., 1.28 times the sound speed), and the reflection height 
A = IO//0. This being a model for a Be star, the azimuthal 
mode is assumed to be retrograde, with i = m = 2 (Rivinius et 
al. 2003), and thus k,Ho = -0.0035. 
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Fig. 5. — Radial dependence of atmospheric parameters for example B2 V 
star: (a) mean density pg modeled with wave pressure (solid line) and without 
(dotted line), (b) bulk wind velocity wq (dashed line), radial wave velocity 
amplitude wi (solid line), sound speed c.i (dot-dashed line), and ratio of res- 
onant wave frequency to local acoustic cutoff frequency (dotted line), (c) 
bulk azimuthal speed uq modeled with wave pressure (solid line) and without 
(dotted line), and the local Keplerian rotation speed (dashed line). 

Figure 5a shows how the inclusion of wave pressure in- 
creases the scale height and gives rise to a significantly shal- 
lower density gradient. The velocity amplitude wi of the res- 
onantly excited wave is shown in Figure 5b, and the initial 
quasi-exponential rise (eq. ifTTl ) is halted by shock dissipa- 
tion. The radial dependence of the time-steady wind speed wo 
is essentially the reciprocal of the density. Figure 5b shows 
that for this model wq just begins to exceed the sound speed 
Cj at the top of the spatial grid, but it is also evident that with- 
out wave pressure this would have occurred at a much lower 
height (z/Hq ^ 10). The angular momentum deposition from 
gv/p.4, is evident in Figure 5c, which shows the abrupt appear- 
ance of nonzero wave pressure beginning at the height where 
shocks first occur. The angular momentum transport satu- 
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rates as uq approaches V^ep because in that limit geff and 

The critical rotation speed Vait is not the same quantity as 
the Keplerian speed Vkcp- The latter is the rotation speed 
needed for a given star to be able to efficiently spin material 
up into orbit and form a disk. Formally, both Vait and Vkcp 
are defined as the rotation speeds required to balance gravity, 
but Vcrit presupposes that the whole star expands out to a crit- 
ical equatorial radius of /J^. « 1 .5 Rp. In the numerical models 
presented here, R^, is held fixed at the value determined by 
the bulk stellar rotation (at and below the photosphere) and 
Vkep = igry^^ is the speed required to balance gravity at that 
radius. For this star, the condition Vrot/Vciit = 0.70 corresponds 
to a ratio Vi-ot/VKep = 0.665. At the photosphere, mq = 335 km 
s~' and Vkbp = 504 km s~', and thus the atmosphere must spin 
up by at least lOc, in order to form a Keplerian disk. Figure 5 
shows that this model contains sufficient angular momentum 
transport to accomplish this spinup below the sonic point of 
the radial outflow. 

5.2. Varying the Stellar Parameters 

The model described above serves as an illustration of how 
resonant excitation, shock dissipation, and wave pressure may 
act together to spin up the outer layers of a rotating star. This 
outcome, though, depends on several key parameters that we 
vary below to explore whether this process can truly be an 
explanation for observed Be-star disks. 

Figure 6 shows contours of various scalar properties taken 
from a two-dimensional grid of models that vary both the ro- 
tation speed Viot (vertical axis) and the vertical NRP veloc- 
ity amplitude Vpuk (horizontal axis) at the photosphere of the 
B2 V star The grid contained 150 values of Viot/Vcrit dis- 
tributed linearly between 0.01 and 0.995, and 150 values of 
Vpuis distributed logarithmically between 0.1 and 100 km s"'. 
All other parameters were held fixed at the values given in 
§5.1. The quantities shown in Figures 6a-c are measured 
at the top of the spatial grid (z = 80//o)- Figure 6a shows 
whether the wave pressure was able to increase the bulk az- 
imuthal speed mq to the local Keplerian speed Vkep by the top 
of the grid. It is clear that a star requires both a relatively high 
photospheric rotation speed (so that mq does not have so far to 
go) and a substantial pulsation amplitude to deposit sufficient 
angular momentum. However, this appears to be able to oc- 
cur for some cases when Vi-ot is as small as about 0.6Vci-it (for 
larger NRP amplitudes), or when Vpuis is as small as 1 or 2 km 
s"' (for larger rotation rates). Figure 6b shows the radial wind 
speed Wo in units of the adopted terminal speed Woo- There 
is a rough correlation between efficient spinup to Keplerian 
rotation and a low wind speed at the top of the grid. 

Figure 6c shows contours of the mean density po at the top 
of the radial grid, which we take as a proxy for the density 
at the "inner edge" of the decretion disk. (For the models 
exhibiting Keplerian rotation, the density gradient is shallow 
and the plotted values of po are insensitive to the exact height 
assumed for the upper edge of the grid.) The densities in Fig- 
ure 6c range between a minimum of 3 x 10""' g cm"-' (lower 
left) and a maximum of 8 x 10"'' g cm"-' (upper right), with 
the values exceeding 10"''* g cm"^ tending to occur in the 
regime of parameter space filled by Keplerian disks. This 
range is consistent with observed inner disk densities. Tradi- 
tional measurements from infrared excess (Waters et al. 1987) 
and visible linear polarization (McDavid 2001) give values 
between about 10"'^ and 1.5 x 10"" g cm"'. More recent 
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Fig. 6. — Properties at the largest height (z = 8O//0) for a grid of B2 V 
star models that vary Vrot and Vpui, independently of one another: (a) ratio 
of rotation speed uq to Keplerian orbital speed VKep, Q>) ratio of radial wind 
speed n'o to assumed terminal speed iVoc, and (c) mean mass density po. 
Labels for constant contour values are given next to each curve. Panel (d) 
shows the heuristic dividing line between winds and disks (see text). See the 
last page of this paper for the full two-column version of this figure. 

determinations that make use of Ha profiles and interferom- 
etry (e.g., Gies et al. 2007; Jones et al. 2008) exhibit values 
that vary a bit more widely; essentially from ~10"'^ to 10"^ g 
cm"^. The latter end of this range, however, overlaps with es- 
timated photospheric densities, and thus these values may be 
inconsistent with observational evidence that the photospheric 
layers rotate more slowly than Vkcp- 

Figure 6d illustrates the "bifurcation" of the models be- 
tween winds and disks. For the purposes of this dividing line, 
a model is considered to have produced a Keplerian disk only 
if both Mo/^Kep > 0.95 and wq/w 00 < 0.1. The latter condi- 
tion ensures that an equatorial wind has not yet accelerated 
significantly by the top of the grid, which would invalidate 
the neglect of supersonic terms in equation ( l30b . 

Figure 7 shows a similar set of contours as Figure 6a, but 
the horizontal axis varies the reflection height A as a free 
parameter (see § 3.2). The two-dimensional grid contained 
120 values of Vrot/Vcrit distributed linearly between 0.01 and 
0.995, and 200 values of A/Z/q distributed logarithmically be- 
tween 1 and 500. This figure also includes the bifurcation 
curve as defined for Figure 6d. The photospheric NRP veloc- 
ity amplitude was held fixed at 20 km s"'. The combined 
effects of resonant excitation, shock dissipation, and wave 
pressure appear to produce roughly the same amount of an- 
gular momentum transport for 1 < A/Z/q < 20. The "kink" 
in the contours at A/i/o ~ 3.5 corresponds to the reflection 
height at which the efficiency of resonant excitation saturates 
to unity somewhere above the photosphere (see Fig. 3). The 
largest values of the reflection height (A/ Ho > 60) do not give 
enough resonant excitation to produce substantial angular mo- 
mentum transport in the upper atmosphere, so the contours 
approach the same kind of limiting shapes as are seen in Fig- 
ure 6 for the limit Vpuis ^ 0. 

The final grid of models was produced by varying the rota- 
tion rate and the spectral type of the star, while keeping both 
Vpuis and A fixed at their standard values of 20 km s"' and 
IO//0, respectively. The grid contained 200 values of Vrot/Vcrit 
distributed linearly between 0.01 and 0.995, and 55 values of 
the spectral type in half-subtype increments between 03 and 
FO. Figure 8 shows contours of both mq/Vkcp (Fig. 8fl) and 
wq/woo (Fig. 8/7) along with the wind/disk bifurcation curve 
that is constrained by the combination of these two quantities. 
As before, the stellar parameters for each spectral type were 
taken from the and Teff calibration of de Jager & Nieuwen- 
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Fig. 7. — Contours of no/Vkep plotted similarly to Figure 6a, but for a grid 
of B2 V star models that vary Vi-ot and the reflection height A independently 
of one another. The photospheric pulsation amplitude is held fixed at 20 km 
s"' . The disk/wind boundary from Figure 6d is also shown {dotted line). 

huijzen (1987), and M* was determined from the evolutionary 
tracks of Claret (2004). Fits for the mass-luminosity relations 
for luminosity classes III, IV, and V were given by Cranmer 
(2005). The models shown in Figure 8 follow the main se- 
quence (i.e., luminosity class V). 

The position of the dividing line between winds and disks 
depends (somewhat) on assumptions made about the equato- 
rial mass loss rate and terminal speed. The fitting relations of 
Vink et al. (2000) were used to estimate M as a function of the 
stellar parameters. The terminal speed was assumed to scale 
as Woo = l-8Vesc, with the escape speed Vesc calculated at the 
pole (i.e., to estimate the polar wind's terminal speed). The 
choice of the numerical constant 1.8 was motivated to repro- 
duce the value of w^o = 1500 km s"' used in the standard B2 V 
model discussed above. Although this constant has been seen 
to be as large as ^3 for earlier O stars (Prinja et al. 1990), the 
shapes of the contours in Figure 8 were not seen to change 
when this constant was varied between 1.8 and 3. 

Figure 8 pinpoints early B-type stars as those most likely to 
form Keplerian decretion disks with largest range of possible 
photospheric rotation rates. The O-type stars have stronger 
winds that require a more self-consistent (subsonic to super- 
sonic) treatment of the radial momentum equation. Late B- 
type and all A-type stars appear to require only slightly larger 
rotation rates than early B stars to produce disks, but one 
should recall that the NRP amplitude was held fixed at 20 km 
s"' across these spectral types. In fact, most "normal" A stars 
(along and above the main sequence) are not as pulsationally 
unstable as the j3 Cep and SPB stars, so their NRP amplitudes 
are likely to be significantly smaller than assumed here. 

6. DISCUSSION 

6.1. Time Steady Disk Formation 

How believable are the models presented above? This is a 
nontrivial question to ask when attempting to solve a problem 
that has plagued a discipline for three quarters of a century. 
The chain of events described in this paper remains some- 
what speculative mainly because of the need to assume a fi- 
nite reflection time (or reflection height A) that provides a 
time-steady yield of resonantly excited waves. If the resonant 
waves are much weaker than those predicted with A < 2QHq 
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to Figures 6a~h, but for a grid of main sequence models that vary Viot and the 
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(see Fig. 7), then this pulsational mechanism may not work. 
Additionally, if the frequency of resonant oscillations varies 
with height to stay equal to the local acoustic cutoff frequency 
(rather than remain at the piston's cutoff frequency), then the 
waves may not propagate upward. In that case, the phase shift 
factors $1 and $2 (eq. 1271 ) could be much smaller than if the 
waves were propagating, and the ability to transport angular 
momentum would be reduced. 

It is important to note, though, that there were several places 
in the above analysis where the wave coupling effects were 
decidedly underestimated.'^ First, equation ( fTTT i neglects all 
higher-order transient terms (i.e., t~^l^, t~''l'^, and so on) in the 
series expansion for the resonantly excited wave power Sec- 
ond, equation ( |23] | uses the local scale height Hiz) in the de- 
nominator of the wave pressure, rather than the smaller value 
at the piston that would have been more consistent with the 
use of a constant frequency uj. Third, the Eulerian mass con- 
servation equation (eq. 1281 ) was solved for wq instead of the 
Lagrangian version, resulting in a potentially weaker impact 
on the angular momentum transport. Fourth, this analysis to- 

In part, this was done to counteract some of the "irrational exuberance" 
that an advocate for a new theoretical model is likely to exhibit! 
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tally neglected the fact that even purely evanescent waves can 
propagate energy in the horizontal direction (Walterscheid & 
Hecht 2003), and thus they may be able to exert a substantial 
wave pressure force in the azimuthal direction — even if they 
do not convert any of their energy into the resonant mode. 

This paper made the implicit assumption that many early- 
type Be stars are rotating significantly below the critical limit, 
as suggested by, e.g., Slettebak (1982), Mennickent et al. 
(1994), Porter (1996), and Cranmer (2005). If instead all Be 
stars are rotating close to critically (i.e., Vrot/Krit ^ 0.95) as 
suggested by Maeder & Meynet (2000) and Townsend et al. 
(2004), the formation of disks could be accomplished with 
much more inefficient — and less speculative — physical mech- 
anisms (see also Owocki 2005). In any case, for Be stars un- 
dergoing NRP with amplitudes of order the sound speed in 
their photospheres, many of the processes described in this 
paper should be occurring nevertheless. 

6.2. Long-Term Be/Bn Phase Variability 

The most prominent feature of Be stars that has not yet been 
addressed here is their substantial variability. Most notice- 
ably, Be stars are known to exhibit phase changes between 
normal B (i.e., "Bn"), standard Be, and Be-sheU properties 
(e.g., Slettebak 1988; Kogure 1990). These changes are typi- 
cally interpreted as variations in the density and geometry of 
the circumstellar disk, though long-term changes in the prop- 
erties of the photosphere have also been inferred (Doazan et 
al. 1986). 

Although it was suggested earlier (§ 1) that a collection 
of other proposed mechanisms for producing Be-star disks 
may be responsible for some of the observed variability, it 
is worthwhile to also examine whether NRP/wave coupling 
would also contribute. In the models presented in this paper, 
the support of Be-star disks is provided ultimately by the ki- 
netic and thermal energy in NRP motions. Thus, variations in 
the overall properties of the disk could be caused by variations 
of the stellar pulsations on the same kinds of time scales (i.e., 
long-term changes in the NRP "piston"). In fact, rapid stellar 
rotation has been shown to be coimected with several kinds of 
pulsational variability and stochasticity: 

1. The rotational splitting of NRP modes gives rise to 
closely separated periods that may be responsible for 
beating and amplitude modulations on time scales of 
weeks to months. For example, Rivinius et al. (1998) 
found a set of periods for /x Cen that are separated 
by fractional frequency differences of order 0.01-0.02, 
such that the beating maxima and minima can occur 
with periods up to several years. There is also long-term 
drift in the precise values of many observed frequencies 
that could indicate beating or other kinds of multi-mode 
interactions (e.g., Stefi & Balona 1996; Sareyan et al. 
1998). It has also been suggested (Henrichs 1984) that 
rotational splitting — especially for rapid rotation — ^may 
create additional "channels" of energy transfer between 
NRP and the mean state, and may even allow more pul- 
sational energy to escape from the star. 

2. As NRP modes grow to large amplitudes (i.e., Vpuis ^ 
Cs) a host of nonlinear wave-wave interactions be- 
come possible that are not present in the linear, small- 
amplitude limit. There is observational evidence for a 
given NRP mode to grow slowly in amplitude, seem- 
ingly at the expense of another mode (Smith 1986). As 



a pulsation mode becomes nonlinear, it may even lose 
sinusoidal coherence by "interfering with itself" as it 
circumnavigates the star in longitude (e.g.. Smith et al. 
1987). If this occurs it could lead to shock steepening 
and dissipation in (j), instead of in r as assumed in § 3.3. 
The passage of an azimuthally propagating shock past a 
given point on the star could also act as a nonlinear trig- 
ger to restart the wave reflection "clock" for resonant 
excitation. Osaki (1999) suggested that a combination 
of these effects — along with feedback from the dense 
circumstellar envelope — ^would also lead to the even- 
tual quenching of angular momentum transport, and 
thus to the dissipation of the disk. Such a repeated 
relaxation-oscillation cycle could be responsible for the 
phase variations between no disks (Bn), standard disks 
(Be), and optically thick disks (Be-sheU). 

3. The existence of subsurface convection could both im- 
pact the long-term growth or decay of pulsation modes 
and act as a seed for circumstellar variability (e.g.. Cox 
1980; Balmforth 1992; Gabriel 1998; Cantiello et al. 
2009). Traditionally, B-type stars are not believed to 
have subsurface convection zones. There is some indi- 
rect evidence, however, that strong turbulent motions 
may be present near the surfaces of rapidly rotating 
early-type stars (Smith 1970; Kodaira 1980; Dolginov 
& Urpin 1983). Some recent models of oblate rapid 
rotators do show the emergence of a near- surface con- 
vection cell at the equator (Espinosa Lara & Rieutord 
2007; Maeder et al. 2008). Even nonrotating hot stars — 
sufficiently above the main sequence — may exhibit en- 
hanced subsurface convection due to newly discovered 
opacity peaks (Cantiello et al. 2009). 

Although strictly unrelated to the stellar pulsations, an- 
other major source of long-term Be-star variability appears 
to be the naturally slow precession of non-axisymmetric den- 
sity fluctuations through the dense Keplerian disk (e.g., Pa- 
paloizou et al. 1992; Savonije & Heemskerk 1993; Fift & 
Harmanec 2006). Gas parcels rotating under the influence of 
a non-point-mass gravitational potential undergo elliptical or- 
bits with slow (decade timescale) pattern speeds. The modes 
that grow the fastest seem to be m = 1, or "one-armed," in- 
stabilities that appear to match the observed properties of the 
evolving violet-to-red (VIR) ratio of double-peaked B aimer 
emission line profiles. The ability to produce accurate and 
predictive models of these variations depends on understand- 
ing the origin and supply of angular momentum at the inner 
boundaries of the disks. 

Finally, magnetic fields should not be neglected as an addi- 
tional source of variability. Recent observational studies are 
beginning to reveal a distinct subclass of magnetic Be stars, 
which includes 7 Cas and at least a half-dozen others around 
spectral types BOe-Ble (Smith et al. 2004; Smith & Balona 
2006; Rakowski et al. 2006; Lopes de Oliveira et al. 2007). 
Magnetic interactions between the star and the disk could be 
the source of dynamo action; i.e., another potential source of 
year- to decade-long variations in the circumstellar emission. 
Note that these processes may be acting whether or not the 
stellar wind is channeled along magnetic flux tubes to feed 
gas into the disk (e.g., Cassinelli et al. 2002; Brown et al. 
2004, 2008; Ud-Doula et al. 2008). 
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7. CONCLUSIONS AND FUTURE PROSPECTS 

The primary aim of this paper has been to explore and test a 
set of physical processes that may be responsible for the pro- 
duction of Keplerian decretion disks around classical Be stars. 
Although nonradial pulsations correspond to low-frequency 
evanescent waves in the photospheres of hot stars, they have 
been suggested to give rise to higher frequency resonant os- 
cillations at the acoustic cutoff frequency. If these resonant 
modes grow in amplitude with increasing height, begin to 
propagate upwards, and steepen into shocks, the resulting dis- 
sipation would create substantial wave pressure that both in- 
creases the atmospheric scale height and transports angular 
momentum upwards. Using a reasonable assumption for the 
efficiency of resonant wave generation, this chain of events 
was found to be able to create the inner boundary conditions 
required for a dense Keplerian disk, even when the underly- 
ing photosphere is rotating as slowly as ~60% of its critical 
rotation speed. 

It is important to note that in the models presented above, 
the "spinup" of the atmosphere occurs high enough above the 
photosphere so that standard measurements of the stellar rota- 
tion speed (i.e., vsin; from absorption line widths) would not 
be sensitive to these motions. The density at which the spinup 
is complete occurs several orders of magnitude below the pho- 
tospheric (t « 1) density, and these values are consistent with 
inner-edge disk densities as measured by infrared excess and 
linear polarization techniques (see § 5.2). One way of po- 
tentially verifying that spinup is occurring in this thin layer 
would be to improve on measurements of spectral lines that 
are formed at a range of heights above the photosphere (see, 
e.g., Chen et al. 1989). In addition, the proposed scenario for 
resonant wave excitation implies that the dominant frequency 
of variabiUty should shift from low (evanescent piston) to high 
(acoustic cutoff) values as the height increases. 

In order to further test and verify that NRP-induced wave 
pressure is responsible for spinning up Be-star disk material, 
the physics of the resonant excitation process needs to be sim- 
ulated more accurately. It may be possible to begin with one- 
dimensional simulations (similar to those of Fleck & Schmitz 
1991) to explore both the time-steady yield of energy con- 
version from evanescent to resonant waves and the gradual 
increase in scale height due to radial wave pressure. Also, an 
analytic non-WKB treatment of partial wave reflection may 
shed light on how much intrinsic "leakage" the evanescent 
waves may undergo, even in the absence of resonant conver- 
sion (see Heinemann & Olbert 1980; Cranmer & van Balle- 
gooijen 2005; Cranmer 2008). 

Testing the full set of processes (e.g., resonant excitation, 
wave pressure, and angular momentum transfer) appears to 
require at least two-dimensional simulations of the fluid prop- 
erties in the equatorial plane of a rotating hot star. These 



may take the form of next-generation extensions of studies by 
Kroll & Hanuschik (1997) and Owocki & Cranmer (2002), 
who showed that a large-amplitude NRP-Uke lower bound- 
ary could give rise to intermittent angular momentum trans- 
fer. Simulations may also benefit from an improved treatment 
of viscosity, which would be needed to extend the models up- 
ward to larger distances in the disk midplane (see, e.g., Lee 
et al. 1991; Okazaki 2001). It is likely that once simulations 
reveal the key physical processes, semi-analytic models like 
those presented in this paper can be improved to be more ro- 
bust, predictively accurate, and less dependent on free param- 
eters. 

For the stellar models in which a supersonic wind begins to 
accelerate, a more physically realistic treatment of the radia- 
tive acceleration is needed. For example, the artificial impo- 
sition of a specific value for M should be replaced with a self- 
consistent calculation of the disk's mass loss rate, which may 
be completely uncoupled from that of the polar wind. The 
proper modeling of B-type stellar outflows may also require 
an explicit treatment of collisional coupling between the var- 
ious atomic and ionic species (Babel 1995; Krticka & Kubat 
2001; Owocki & Puis 2002; Votruba et al. 2007) and a consis- 
tent solution for the internal energy equation, including both 
shock heating (e.g., Struck et al. 2004) and radiative scattering 
effects (Gayley & Owocki 1994). 

Lastly, it should be emphasized that future work must in- 
volve not only increased physical realism for the models, but 
also quantitative comparisons with observations. It would 
be beneficial to apply the methodology outlined in this pa- 
per to a set of well-observed stars, rather than to the ide- 
alized spectral- type sequence illustrated in Figure 8. Mea- 
sured stellar properties — including the dominant NRP veloc- 
ity amplitudes — would be used as input constraints for models 
sinoilar to those described in this paper. These models could 
then produce specific predictions for whether there should (or 
should not) be an equatorial disk for each star. Other mod- 
els of the Be phenomenon, such as those summarized in § 1, 
will likely result in a different division of predicted disks ver- 
sus non-disks for the same database of stars. Comparisons of 
each set of predictions with the observed occurrence of disks 
would act as a clear test of the viability of the proposed phys- 
ical processes. 
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APPENDIX 

DERIVATION OF SECOND ORDER WAVE PRESSURE TERMS 

For simplicity, the equations of hydrodynamic mass and momentum conservation are written below in Cartesian, plane-parallel 
coordinates. The vertical direction is z and the relevant horizontal direction for oblique wave propagation is x. This geometry 
is applied to the equatorial plane of a rotating system, where z points radially outward and x points in the prograde azimuthal 
direction (i.e., increasing x means increasing (f) in the direction of rotational motion). We assume symmetry above and below the 
equatorial plane, so that all variations in colatitude (or Cartesian coordinate y) are ignored. 

The two nontrivial velocity components are defined as m = V;^ and w = v^. The equations of mass and momentum conservation 
are written as 

^ + ^{pu) + ^{pw) = Q (Al) 
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d d d dP 

— (pw) + — (puw) + —(pv/) + — -pg = Q ( A2) 

at ox oz oz 

d d J d dP 

It is important to write the above equations in their full "conservation form." Let us assume that all dynamical variables can be 
split into the sum of zeroth order (time steady) and first order (sinusoidally oscillating) components. Initially, the radial wind 
velocity wo(z) and the mean rotation velocity uq(z) are assumed to have arbitrary dependences on height. The mean density scale 
height is defined as H = -po/idpo/dz). All zeroth order properties are assumed to be constant in both time f and horizontal 
position X. 

Retaining only the zeroth order terms yields equations for the mean atmospheric stratification in the absence of wave pressure. 
Collecting the first order terms yields the dispersion relation as well as relative amplitude relations for linear waves. These 
are discussed by, e.g., Mihalas & Mihalas (1984) and Wang et al. (1995), and the first-order expressions themselves are given 
in equations (7.82)-(7.90) of Cranmer (1996). The remainder of this Appendix deals with the second order terms; i.e., those 
that involve products of two first order fluid quantities. If the time average of a given second order term (taken over either a 
single wave period of an integer multiple of periods) does not vanish, then it gives rise to a net contribution to the time-steady 
conservation equations. 

For example, the period-averaged mass flux conservation equation can be written as 

^(P()W())+^^(piMi)^ + (^^(piWi)^ =0. (A4) 

The period-averaging is denoted by angle brackets. These terms can be simplified further by exploring some of the mathematical 
properties of such time averages. Consider a complex first order quantity / = (/,. + ///)e'*"'"*' '\ where the real amplitudes /,. and /,■ 
are assumed to contain an unspecified dependence on z- For two such quantities (/ and g), the Eulerian average of their product 
over a period Itt/lu has the following real part: 

ifg) = \ifr8r + f:gd = \(f*g + fg*) (A5) 

where /* denotes the complex conjugate of / and this averaged quantity is constant in t and x. An additional and useful identity 
can be derived for these kinds of sinusoidal time variations, 

and a similar identity can be defined by replacing f with x, since the x variations are ideally sinusoidal as well. (This cannot be 
done for the variations in z, however.) Thus, the second term in equation (IA4b is zero, and the radial derivative in the third term 
can be taken outside the angle brackets to obtain 

d 

— {p()Wo+ {p\wi)) = . (A7) 

This is the Lagrangian form of the mass conservation equation, which indicates that longitudinal and compressive waves can be 
responsible for a nonzero mass flux even in the absence of a mean local (Eulerian) flow velocity wq. The ramifications of the 
(piWi) term (also sometimes called the "Stokes drift") are discussed in § 4. 

The sum of the second order terms in the x and z components of the momentum conservation equation form a vector that we 
define as S = Sx^x + S^e^. Dimensionally, S has the units of a pressure gradient, or force per unit volume, and its components are 

■^-^ = ^(PlMl)+^(POM? + 2/OlMoMl)+^(POMlWl+PlMoWl+PlMlWo) (AS) 

ot ox Oz 

d d d 

These expressions can be simplified by noticing that the derivatives with respect to t and x can be neglected when constructing 
the period averages of Sx and S^. Thus, the mean azimuthal wave pressure gradient is given by 

{Sx) = (^^ipouiwi)^ + (^-^(uopiwi)^ + (^-^(wopiui)^ . (AlO) 

The first and second terms in equation ( lAlOl ) are referred to by Lee & Saio (1993) and Lee (2007) as occurring on transport time 
scales Ti and T2, respectively. The third term depends on the Eulerian mean velocity wo and is traditionally ignored in stellar 
interiors models (for which the only mean motions are meridional circulations with speeds many orders of magnitude smaller 
than the pulsation amplitudes). We also neglect the wq term here (see below), and the remaining terms can be written as 

{Sx) = ^ (po(miVI'i)) + ^ (mo(piWi)) . (All) 

Because the first term above depends on ui, its overall sign thus depends on the sign of kx (see eq. IQ). The second term, though, 
has the same sign no matter whether the waves are prograde or retrograde. For this term, it is straightforward to see that when 
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Pi and wi are in phase with one another (e.g., for an upwardly propagating wave) and uq is positive (i.e., for bulk rotation in the 
+x direction), it behaves similarly to the radial wave pressure. In other words, if the pressure-like quantity being differentiated 
with respect to z decreases as one goes up in height, it will result in a positive force (i.e., an upward acceleration for and an 
acceleration in the prograde direction for S^). 
The vertical (radial) wave pressure gradient can be condensed down to only two terms, with 



(- 

\dz 



+ 2 



(wopiwi) 



(A12) 



(A13) 



As above, the term depending explicitly on wo will be neglected here. Thus, the single remaining term 

(5e) = ^(po(w?)) 

is equivalent to the plane-parallel limit of the wave pressure gradient derived elsewhere (e.g., Jacques 1977; Koninx 1992) in 
the limit of 7 = 1 (isothermal fluctuations). Also, the above expression is equivalent to equation ( l22b as long as the waves are 

sinusoidal in shape (s = 2). 

It should be emphasized that the above form of the radial wave pressure (eq. IIA13I ) is the same as that deriv ed by Jacques 
(1977) and others, despite the fact that they do not make the assumption that wo = 0. In fact, they apply equation (IA13b directly 
to the problem of how waves can help accelerate supersonic winds. Thus, it is likely that taking the limit of wo ^ in the above 
equations is equivalent to transforming into the local comoving frame of reference of the fluid. ^ The proper Reynolds stresses in 
Sx and are those that add momentum to the fluid in the comoving (Lagrangian) frame. 
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Fig. 6. — Properties at the largest height (z = SOHq) for a grid of B2 V star models that vary Vrot and Vpuis independently of one another: (a) ratio of rotation 
speed Mo to Keplerian orbital speed Vkep, (b) ratio of radial wind speed wq to assumed terminal speed Woo, and (c) mean mass density po- Labels for constant 
contour values are given next to each curve. Panel (d) shows the heuristic dividing line between winds and disks (see text). 



